### REPLICATION FILE -- APPENDIX
### Jonathan Homola
### "The Effects of Women's Descriptive Representation on Government Behavior"
### Legislative Studies Quarterly

## clear environment, set seed, install/load packages
rm(list=ls())
set.seed(12435); options(stringsAsFactors=F)
# install.packages("foreign"); install.packages("stargazer")
# install.packages("survey"); install.packages("xtable")
# install.packages("data.table"); install.packages("pBrackets")
library(foreign); library(stargazer); library(survey); library(psych); library(xtable); library(data.table); library(pBrackets)

# GLM confidence intervals
glm.cis <- function(preds, ses, alpha, df){
  tval = qt((1-alpha)/2, df, lower.tail = F)
  raw_conf = cbind(preds-(tval*ses), preds+(tval*ses))
  trans_conf = exp(raw_conf)/(1+exp(raw_conf))
  out <- cbind(preds, raw_conf, exp(preds)/(1+exp(preds)), trans_conf)
  colnames(out) <- c('eta', 'low.eta', 'hi.eta', 'pred', 'low.pred', 'hi.pred')
  return(out)
}

# calculate R^2 for survey design models
fit.svyglm = function(svyglm, digits=3){
  if(methods::is(svyglm, "svyglm")!=TRUE) message("Warning: Not a svyglm object.")
  r2 = (svyglm$null.deviance - svyglm$deviance) / svyglm$null.deviance 
  adjust = svyglm$df.null / svyglm$df.residual
  value  = 1 - ((1-r2) * adjust)
  results = c(R2=round(r2, digits), adjR2=round(value, digits))
  names(results) = c("R-Squared","     Adjusted R-Squared")
  return(results)
}

## set working directory
setwd(" ... ")

## read in the dataset 
pledges_data <- read.dta("Homola_Pledges.dta")

# create subset variable for Table 2 in Thomson et al
pledges_data$subset1 <- ifelse(pledges_data$govparty==1 & pledges_data$sq==0 & (pledges_data$excllink==0 | is.na(pledges_data$excllink)==T), 1, 0)
df_tab2 <- subset(pledges_data, subset1==1)

# survey weights and design for full data set
design.cab <- svydesign(ids=~idmanifesto, weights=~sampwtmod7770, data=df_tab2)


###
### Table A.2: Descriptive statistics
###

xtable(describe(df_tab2[,c("femaleleader_complete", "women_cabinet",
                        "fulfil2", "gtsinmin", "gtcomaj", "gtcomin", "chex", "ministry", "govlogrilerange",
                        "HHgovpr", "presidential", "semipres", "bicameral", "federal", "EUmember", "growav",
                        "xtimeyrs", "opppreel", "outalways", "npledgediv10", "precoagree", "degomedlogrile",
                        "y1980s", "y1990s", "y2000s", "subset")])[,c(2,3,4,8,9)])


## models for Table 1 / Appendix Table A.3
m1a <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m1b <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m2a <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m2b <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m3a <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)


###
### Table A.3: Women's descriptive representation and pledge fulfillment -- full results
###

stargazer(m1a, m2a, m3a, m1b, m2b, m3b,
          covariate.labels=c("Single-party minority", "Coalition majority", "Coalition minority",
                             "Chief executive", "Relevant portfolio", "Ideological range",
                             "Herfindahl index", "Presidentialism", "Semi-presidentialism",
                             "Bicameralism", "Federalism", "EU member", "GDP growth",
                             "Duration in years", "Opposition parties with experience",
                             "Opposition parties without experience", "Number of pledges (/10)",
                             "Pre-election coalition", "Ideological distance to median legislator",
                             "1980s", "1990s", "2000s", "Subset of pledges tested",
                             "Female party leader", "Share of women in cabinet", "Constant"),
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor",
          add.lines = list(c("Thomson et al. covariates", "checkmark", "checkmark", "checkmark", "checkmark", "checkmark", "checkmark"),
                           c("Party family FEs", "", "checkmark", "", "", "checkmark", ""),
                           c("Party FEs", "", "", "checkmark", "", "", "checkmark")))

## check r2
fit.svyglm(m1a)
fit.svyglm(m2a)
fit.svyglm(m3a)
fit.svyglm(m1b)
fit.svyglm(m2b)
fit.svyglm(m3b)

## check number of party programs
aggregate(df_tab2[, 5], list(df_tab2$idmanifesto, df_tab2$year), mean)

## check variation on female party leader with partyid FEs
part1 <- aggregate(df_tab2[, 83], list(df_tab2$partyid, df_tab2$year), mean)
aggregate(part1[, 3], list(part1$Group.1), mean)
## only two parties with change in femaleleader: UK Conservatives and Irish Progressive Democrats


###
### Figure A.1: Predicted probabilities of pledge fulfillment for different levels of the share of women in cabinet -- party fixed effects model
###

# wcab, party FE (partyid=702, german CDU/CSU)
# drop covariates that are estimated as NA because of collinearity: presidential + semipres + bicameral + federal 
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
newdat_wcabiq <- data.frame(gtsinmin = rep(0, 2), gtcomaj = rep(0, 2), gtcomin = rep(0, 2),
                            chex = rep(1, 2), ministry = rep(1, 2), 
                            govlogrilerange = rep(0.4061, 2), HHgovpr = rep(0.7503, 2),
                            presidential = rep(0, 2), semipres = rep(0, 2), bicameral = rep(1, 2),
                            federal = rep(0, 2), EUmember = rep(1, 2), growav = rep(2.455, 2),
                            xtimeyrs = rep(3.707, 2), opppreel = rep(0, 2), outalways = rep(0, 2),
                            npledgediv10 = rep(17.74, 2), precoagree = rep(0, 2), 
                            degomedlogrile = rep(0.2491240, 2), y1980s = rep(0, 2),
                            y1990s = rep(0, 2), y2000s = rep(1, 2), subset = rep(0, 2),
                            women_cabinet = as.numeric(summary(df_tab2$women_cabinet)[c(2,5)]),
                            partyid = rep(factor(702, levels(as.factor(df_tab2$partyid))), 2))
newdat_wcabf <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                           chex = rep(1, 200), ministry = rep(1, 200), 
                           govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                           presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                           federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                           xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                           npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                           degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                           y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                           women_cabinet = seq(0, 0.6, length.out=200),
                           partyid = rep(factor(702, levels(as.factor(df_tab2$partyid))), 200))
preds_wcabiq <- predict(m3b, newdata = newdat_wcabiq, se.fit = TRUE)
preds_wcabf <- predict(m3b, newdata = newdat_wcabf, se.fit = TRUE)
cis_wcabiq <- round(glm.cis(as.numeric(preds_wcabiq), sqrt(attr(preds_wcabiq, "var")), 0.95, m3b$df.residual),4)
cis_wcabf <- round(glm.cis(as.numeric(preds_wcabf), sqrt(attr(preds_wcabf, "var")), 0.95, m3b$df.residual),4)
iqrange <- cbind(summary(df_tab2$women_cabinet)[c(2,5)],cis_wcabiq[2,4] - cis_wcabf[1,4])

mygray = rgb(153, 153, 153, alpha = 200, maxColorValue = 255)
pdf("cabinet_partyfe.pdf", height=8,width=12)
plot(newdat_wcabf$women_cabinet, cis_wcabf[,4], type="n",xlab="",ylab="",  yaxt="n", xaxt="n",
     ylim = c(0.2,1))
polygon(x = c(newdat_wcabf$women_cabinet, rev(newdat_wcabf$women_cabinet)),
        y = c(cis_wcabf[,5], rev(cis_wcabf[,6])), col = mygray, border = NA)
lines(newdat_wcabf$women_cabinet, cis_wcabf[,4], lwd=2)
rug(df_tab2$women_cabinet)
segments(iqrange[,1], cis_wcabiq[,4], iqrange[,1], rep(0.5,2), lty=2)
segments(iqrange[1,1], 0.5, iqrange[2,1], 0.5, lty = 2)
brackets(iqrange[1,1], 0.4, iqrange[2,1], 0.4, h = -0.05,  ticks = 0.5, lwd=2)
text(abs((iqrange[2,1]-iqrange[1,1])/2)+iqrange[1,1], 0.45, 'Interquartile range', cex=1.5)
text(iqrange[,1], cis_wcabiq[,4]+0.025, round(iqrange[,1],3), cex=1.5)
text(abs((iqrange[2,1]-iqrange[1,1])/2)+iqrange[1,1], 0.3,
     labels=bquote(hat(y)['Q3']-hat(y)['Q1'] ~ '=' ~ .(iqrange[1,2])), cex=1.5)
axis(1, tck=0.03, cex.axis=1.5, mgp=c(0.3, 0.3, 0), lty=1, lwd=0, lwd.ticks = 1)
axis(2, tck=0.03, cex.axis=1.5, mgp=c(0.3, 0.3, 0), lty=1, lwd=0, lwd.ticks = 1, las=2)
title(xlab = 'Share of Women in Cabinet',
      ylab="Probability of Pledge Fulfillment",
      line = 2.25, cex.lab=1.5)
dev.off()


###
### Table A.4: Women's descriptive representation and pledge fulfillment -- including status quo pledges
###

# create subset variable for Table 2 in Thomson et al, include SQ
pledges_data$subset2 <- ifelse(pledges_data$govparty==1 & (pledges_data$excllink==0 | is.na(pledges_data$excllink)==T), 1, 0)
df_tab2b <- subset(pledges_data, subset2==1)
# calculate new weights 
dt <- data.table(df_tab2b)
dt[ , count := .N, by = list(countryn)]
df_tab2b <- as.data.frame(dt)
df_tab2b$weight <- (nrow(df_tab2b)/12)/df_tab2b$count
design.m2b <- svydesign(ids=~idmanifesto, weights=~weight, data=df_tab2b)

# run models
m1a <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab2b, design=design.m2b)
m1b <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2b, design=design.m2b)
m2a <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab2b, design=design.m2b)
m2b <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2b, design=design.m2b)
m3a <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab2b, design=design.m2b)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2b, design=design.m2b)

# create table
stargazer(m1a, m2a, m3a, m1b, m2b, m3b,
          covariate.labels=c("Single-party minority", "Coalition majority", "Coalition minority",
                             "Chief executive", "Relevant portfolio", "Ideological range",
                             "Herfindahl index", "Presidentialism", "Semi-presidentialism",
                             "Bicameralism", "Federalism", "EU member", "GDP growth",
                             "Duration in years", "Opposition parties with experience",
                             "Opposition parties without experience", "Number of pledges (/10)",
                             "Pre-election coalition", "Ideological distance to median legislator",
                             "1980s", "1990s", "2000s", "Subset of pledges tested",
                             "Female party leader", "Share of women in cabinet", "Constant"),
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor",
          add.lines = list(c("Thomson et al. covariates", "checkmark", "checkmark", "checkmark", "checkmark", "checkmark", "checkmark"),
                           c("Party family FEs", "", "checkmark", "", "", "checkmark", ""),
                           c("Party FEs", "", "", "checkmark", "", "", "checkmark")))

## check r2
fit.svyglm(m1a)
fit.svyglm(m2a)
fit.svyglm(m3a)
fit.svyglm(m1b)
fit.svyglm(m2b)
fit.svyglm(m3b)


###
### Table A.5: Women's descriptive representation and pledge fulfillment -- alternative cabinet coding
###

m1c <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinetA, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m2c <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinetA, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m3c <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinetA, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)

stargazer(m1c, m2c, m3c,
          covariate.labels=c("Single-party minority", "Coalition majority", "Coalition minority",
                             "Chief executive", "Relevant portfolio", "Ideological range",
                             "Herfindahl index", "Presidentialism", "Semi-presidentialism",
                             "Bicameralism", "Federalism", "EU member", "GDP growth",
                             "Duration in years", "Opposition parties with experience",
                             "Opposition parties without experience", "Number of pledges (/10)",
                             "Pre-election coalition", "Ideological distance to median legislator",
                             "1980s", "1990s", "2000s", "Subset of pledges tested",
                             "Share of women in cabinet (Alternative)", "Constant"),
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor",
          add.lines = list(c("Thomson et al. covariates", "checkmark", "checkmark", "checkmark"),
                           c("Party family FEs", "", "checkmark", ""),
                           c("Party FEs", "", "", "checkmark")))

fit.svyglm(m1c)
fit.svyglm(m2c)
fit.svyglm(m3c)


###
### Table A.6: Women's descriptive representation and pledge fulfillment -- single-party governments (Thomson et al. [2017], Table 3, Model 1)
###

# create subset variable for Table 3 Model 1
pledges_data$subset3a <- ifelse(pledges_data$coalition==0 & pledges_data$govparty==1 & pledges_data$sq==0 & (pledges_data$excllink==0 | is.na(pledges_data$excllink)==T), 1, 0)
df_tab3a <- subset(pledges_data, subset3a==1)
# weights and design for full data set
design.3a <- svydesign(ids=~idmanifesto, weights=~sampwtmod2946, data=df_tab3a)

# run models
m3a1 <- svyglm(fulfil2 ~ gtminless3yrs + gtminmore3yrs + presidential + semipres + bicameral +federal + growav + opppreel + outalways + npledgediv10 + degomedlogrile + y1980s + y1990s+ y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab3a, design=design.3a)
m3a2 <- svyglm(fulfil2 ~ as.factor(parfam) + gtminless3yrs + gtminmore3yrs + presidential + semipres + bicameral +federal + growav + opppreel + outalways + npledgediv10 + degomedlogrile + y1980s + y1990s+ y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab3a, design=design.3a)
m3a3 <- svyglm(fulfil2 ~ as.factor(partyid) + gtminless3yrs + gtminmore3yrs + presidential + semipres + bicameral +federal + growav + opppreel + outalways + npledgediv10 + degomedlogrile + y1980s + y1990s+ y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab3a, design=design.3a)
m3a4 <- svyglm(fulfil2 ~ gtminless3yrs + gtminmore3yrs + presidential + semipres + bicameral +federal + growav + opppreel + outalways + npledgediv10 + degomedlogrile + y1980s + y1990s+ y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab3a, design=design.3a)
m3a5 <- svyglm(fulfil2 ~ as.factor(parfam) + gtminless3yrs + gtminmore3yrs + presidential + semipres + bicameral +federal + growav + opppreel + outalways + npledgediv10 + degomedlogrile + y1980s + y1990s+ y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab3a, design=design.3a)
m3a6 <- svyglm(fulfil2 ~ as.factor(partyid) + gtminless3yrs + gtminmore3yrs + presidential + semipres + bicameral +federal + growav + opppreel + outalways + npledgediv10 + degomedlogrile + y1980s + y1990s+ y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab3a, design=design.3a)

# create table
stargazer(m3a1, m3a2, m3a3, m3a4, m3a5, m3a6,
          covariate.labels=c("Minority governments of less than 3 years", "Minority governments of at least 3 years", 
                             "Presidentialism", "Semi-presidentialism", "Bicameralism", "Federalism", "GDP growth", 
                             "Opposition parties with experience", "Opposition parties without experience", 
                             "Number of pledges (/10)", "Ideological distance to median legislator", 
                             "1980s", "1990s", "2000s", "Subset of pledges tested", 
                             "Female party leader", "Share of women in cabinet", "Constant"),
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor",
          add.lines = list(c("Party family FEs", "", "checkmark", "", "", "checkmark", ""),
                           c("Party FEs", "", "", "checkmark", "", "", "checkmark")))
fit.svyglm(m3a1)
fit.svyglm(m3a2)
fit.svyglm(m3a3)
fit.svyglm(m3a4)
fit.svyglm(m3a5)
fit.svyglm(m3a6)
## check number of party programs
aggregate(df_tab3a[, 5], list(df_tab3a$idmanifesto, df_tab3a$year), mean)


###
### Table A.7: Women's descriptive representation and pledge fulfillment -- coalition governments (Thomson et al. [2017], Table 3, Model 2)
###

# create subset variable for Table 3 Model 2
pledges_data$subset3b <- ifelse(pledges_data$coalition==1 & pledges_data$govparty==1 & pledges_data$sq==0 & (pledges_data$excllink==0 | is.na(pledges_data$excllink)==T), 1, 0)
df_tab3b <- subset(pledges_data, subset3b==1)
# weights and design for full data set
design.3b <- svydesign(ids=~idmanifesto, weights=~sampwtmod4021, data=subset(df_tab3b, !is.na(sampwtmod4021)))

# run models
m3b1 <- svyglm(fulfil2 ~ gtmajless3yrs + gtminmore3yrs + chex + ministry + govlogrilerange +  HHgovpr + pledgecoagree + bicameral + federal + growav + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab3b, design=design.3b)
m3b2 <- svyglm(fulfil2 ~ as.factor(parfam) + gtmajless3yrs + gtminmore3yrs + chex + ministry + govlogrilerange +  HHgovpr + pledgecoagree + bicameral + federal + growav + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab3b, design=design.3b)
m3b3 <- svyglm(fulfil2 ~ as.factor(partyid) + gtmajless3yrs + gtminmore3yrs + chex + ministry + govlogrilerange +  HHgovpr + pledgecoagree + bicameral + federal + growav + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=df_tab3b, design=design.3b)
m3b4 <- svyglm(fulfil2 ~ gtmajless3yrs + gtminmore3yrs + chex + ministry + govlogrilerange +  HHgovpr + pledgecoagree + bicameral + federal + growav + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab3b, design=design.3b)
m3b5 <- svyglm(fulfil2 ~ as.factor(parfam) + gtmajless3yrs + gtminmore3yrs + chex + ministry + govlogrilerange +  HHgovpr + pledgecoagree + bicameral + federal + growav + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab3b, design=design.3b)
m3b6 <- svyglm(fulfil2 ~ as.factor(partyid) + gtmajless3yrs + gtminmore3yrs + chex + ministry + govlogrilerange +  HHgovpr + pledgecoagree + bicameral + federal + growav + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab3b, design=design.3b)

# create table
stargazer(m3b1, m3b2, m3b3, m3b4, m3b5, m3b6,
          covariate.labels=c("Majority governments of less than three years", "Minority governments of at least 3 years", 
                             "Chief executive", "Relevant portfolio", "Ideological range", "Herfindahl index", 
                             "Agreement between coalition partners", "Bicameralism", "Federalism", "GDP growth", 
                             "Opposition parties with experience", "Opposition parties without experience", 
                             "Number of pledges (/10)", "Pre-election coalition", 
                             "Ideological distance to median legislator", 
                             "1990s", "2000s", "Subset of pledges tested", 
                             "Female party leader", "Share of women in cabinet", "Constant"),
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor",
          add.lines = list(c("Party family FEs", "", "checkmark", "", "", "checkmark", ""),
                           c("Party FEs", "", "", "checkmark", "", "", "checkmark")))
fit.svyglm(m3b1)
fit.svyglm(m3b2)
fit.svyglm(m3b3)
fit.svyglm(m3b4)
fit.svyglm(m3b5)
fit.svyglm(m3b6)
## check number of party programs
df_tab3b_n <- subset(df_tab3b, !is.na(sampwtmod4021))
aggregate(df_tab3b_n[, 5], list(df_tab3b_n$idmanifesto, df_tab3b_n$year), mean)


###
### Table A.8: Women's descriptive representation and pledge fulfillment -- both measures
###

m1a <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m2a <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m3a <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)

stargazer(m1a, m2a, m3a,
          covariate.labels=c("Single-party minority", "Coalition majority", "Coalition minority",
                             "Chief executive", "Relevant portfolio", "Ideological range",
                             "Herfindahl index", "Presidentialism", "Semi-presidentialism",
                             "Bicameralism", "Federalism", "EU member", "GDP growth",
                             "Duration in years", "Opposition parties with experience",
                             "Opposition parties without experience", "Number of pledges (/10)",
                             "Pre-election coalition", "Ideological distance to median legislator",
                             "1980s", "1990s", "2000s", "Subset of pledges tested",
                             "Female party leader", "Share of women in cabinet", "Constant"),
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor",
          add.lines = list(c("Thomson et al. covariates", "checkmark", "checkmark", "checkmark"),
                           c("Party family FEs", "", "checkmark", ""),
                           c("Party FEs", "", "", "checkmark")))

## check r2
fit.svyglm(m1a)
fit.svyglm(m2a)
fit.svyglm(m3a)


###
### Table A.9: Country-election years sorted by Share of women in cabinet
###

desc_table <- aggregate(df_tab2[, 80], list(df_tab2$countryn, df_tab2$year), mean)
desc_table <- desc_table[order(desc_table$x),]
print(xtable(desc_table, digits=c(0,0,0,2)), include.rownames=F)


###
### Table A.10: Country-election years sorted by Pledge fulfillment
###

desc_table <- aggregate(df_tab2[, 10], list(df_tab2$countryn, df_tab2$year), mean)
desc_table <- desc_table[order(desc_table$x),]
print(xtable(desc_table, digits=c(0,0,0,2)), include.rownames=F)


###
### Table A.11: Women's descriptive representation and pledge fulfillment -- without Sweden
###

nosweden <- subset(df_tab2, countryn!="Sweden")
# calculate weights for outlier analysis
dt <- data.table(nosweden)
dt[ , count := .N, by = list(countryn)]
nosweden_weights <- as.data.frame(dt)
nosweden_weights$weight <- (nrow(nosweden_weights)/11)/nosweden_weights$count
design.nosweden <- svydesign(ids=~idmanifesto, weights=~weight, data=nosweden_weights)

m1a <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nosweden, design=design.nosweden)
m1b <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nosweden, design=design.nosweden)
m2a <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nosweden, design=design.nosweden)
m2b <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nosweden, design=design.nosweden)
m3a <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nosweden, design=design.nosweden)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nosweden, design=design.nosweden)
stargazer(m1a, m2a, m3a, m1b, m2b, m3b,
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor")
## check r2
fit.svyglm(m1a)
fit.svyglm(m2a)
fit.svyglm(m3a)
fit.svyglm(m1b)
fit.svyglm(m2b)
fit.svyglm(m3b)
## check number of party programs
aggregate(nosweden[, 5], list(nosweden$idmanifesto, nosweden$year), mean)


###
### Table A.12: Women's descriptive representation and pledge fulfillment -- without the US
###

nous <- subset(df_tab2, countryn!="USA")
# calculate weights for outlier analysis
dt <- data.table(nous)
dt[ , count := .N, by = list(countryn)]
nous_weights <- as.data.frame(dt)
nous_weights$weight <- (nrow(nous_weights)/11)/nous_weights$count
design.nous <- svydesign(ids=~idmanifesto, weights=~weight, data=nous_weights)

m1a <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nous, design=design.nous)
m1b <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nous, design=design.nous)
m2a <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nous, design=design.nous)
m2b <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nous, design=design.nous)
m3a <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nous, design=design.nous)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nous, design=design.nous)
stargazer(m1a, m2a, m3a, m1b, m2b, m3b,
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor")
## check r2
fit.svyglm(m1a)
fit.svyglm(m2a)
fit.svyglm(m3a)
fit.svyglm(m1b)
fit.svyglm(m2b)
fit.svyglm(m3b)
## check number of party programs
aggregate(nous[, 5], list(nous$idmanifesto, nous$year), mean)


###
### Table A.13: Women's descriptive representation and pledge fulfillment -- without the UK
###

nouk <- subset(df_tab2, countryn!="United Kingdom")
# calculate weights for outlier analysis
dt <- data.table(nouk)
dt[ , count := .N, by = list(countryn)]
nouk_weights <- as.data.frame(dt)
nouk_weights$weight <- (nrow(nouk_weights)/11)/nouk_weights$count
design.nouk <- svydesign(ids=~idmanifesto, weights=~weight, data=nouk_weights)

m1a <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nouk, design=design.nouk)
m1b <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nouk, design=design.nouk)
m2a <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nouk, design=design.nouk)
m2b <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nouk, design=design.nouk)
m3a <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + femaleleader_complete, family=quasibinomial(link="logit"), data=nouk, design=design.nouk)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=nouk, design=design.nouk)
stargazer(m1a, m2a, m3a, m1b, m2b, m3b,
          digits=2, star.cutoffs = c(0.05, 0.01, NA), dep.var.labels.include = FALSE, no.space=T,
          dep.var.caption="Outcome variable: Probability of pledge fulfillment",
          omit = "as.factor")
## check r2
fit.svyglm(m1a)
fit.svyglm(m2a)
fit.svyglm(m3a)
fit.svyglm(m1b)
fit.svyglm(m2b)
fit.svyglm(m3b)
## check number of party programs
aggregate(nouk[, 5], list(nouk$idmanifesto, nouk$year), mean)


###
### Figure A.2: Predicted probabilities of pledge fulfillment for different levels of (a) gender of party leader and (b) share of women in cabinet
###

# "linear" effects
m1b <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m2b <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)
m3b <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab)

newdat_wcabf1 <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                            chex = rep(1, 200), ministry = rep(1, 200), 
                            govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                            presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                            federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                            xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                            npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                            degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                            y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                            women_cabinet = seq(0, 0.6, length.out=200))
newdat_wcabf2 <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                            chex = rep(1, 200), ministry = rep(1, 200), 
                            govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                            presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                            federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                            xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                            npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                            degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                            y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                            women_cabinet = seq(0, 0.6, length.out=200),
                            parfam = rep(factor(60, levels(as.factor(df_tab2$parfam))), 200))
newdat_wcabf3 <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                            chex = rep(1, 200), ministry = rep(1, 200), 
                            govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                            presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                            federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                            xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                            npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                            degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                            y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                            women_cabinet = seq(0, 0.6, length.out=200),
                            partyid = rep(factor(701, levels(as.factor(df_tab2$partyid))), 200))
preds_wcabf1 <- predict(m1b, newdata = newdat_wcabf1, se.fit = TRUE)
cis_wcabf1 <- round(glm.cis(as.numeric(preds_wcabf1), sqrt(attr(preds_wcabf1, "var")), 0.95, m1b$df.residual),4)
preds_wcabf2 <- predict(m2b, newdata = newdat_wcabf2, se.fit = TRUE)
cis_wcabf2 <- round(glm.cis(as.numeric(preds_wcabf2), sqrt(attr(preds_wcabf2, "var")), 0.95, m2b$df.residual),4)
preds_wcabf3 <- predict(m3b, newdata = newdat_wcabf3, se.fit = TRUE)
cis_wcabf3 <- round(glm.cis(as.numeric(preds_wcabf3), sqrt(attr(preds_wcabf3, "var")), 0.95, m3b$df.residual),4)


# "non-linear" effects
df_tab2$women_cabinet_sq <- df_tab2$women_cabinet*df_tab2$women_cabinet
design.cab.sq <- svydesign(ids=~idmanifesto, weights=~sampwtmod7770, data=df_tab2)
m1bsq <- svyglm(fulfil2 ~ gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet + women_cabinet_sq, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab.sq)
m2bsq <- svyglm(fulfil2 ~ as.factor(parfam) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + presidential + semipres + bicameral + federal + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet + women_cabinet_sq, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab.sq)
m3bsq <- svyglm(fulfil2 ~ as.factor(partyid) + gtsinmin + gtcomaj + gtcomin + chex + ministry + govlogrilerange +  HHgovpr + EUmember + growav + xtimeyrs + opppreel + outalways + npledgediv10 + precoagree + degomedlogrile + y1980s + y1990s + y2000s + subset + women_cabinet + women_cabinet_sq, family=quasibinomial(link="logit"), data=df_tab2, design=design.cab.sq)

newdat_wcabf1sq <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                              chex = rep(1, 200), ministry = rep(1, 200), 
                              govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                              presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                              federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                              xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                              npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                              degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                              y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                              women_cabinet = seq(0, 0.6, length.out=200),
                              women_cabinet_sq = (seq(0, 0.6, length.out=200))^2)
newdat_wcabf2sq <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                              chex = rep(1, 200), ministry = rep(1, 200), 
                              govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                              presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                              federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                              xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                              npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                              degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                              y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                              women_cabinet = seq(0, 0.6, length.out=200),
                              women_cabinet_sq = (seq(0, 0.6, length.out=200))^2,
                              parfam = rep(factor(60, levels(as.factor(df_tab2$parfam))), 200))
newdat_wcabf3sq <- data.frame(gtsinmin = rep(0, 200), gtcomaj = rep(0, 200), gtcomin = rep(0, 200),
                              chex = rep(1, 200), ministry = rep(1, 200), 
                              govlogrilerange = rep(0.4061, 200), HHgovpr = rep(0.7503, 200),
                              presidential = rep(0, 200), semipres = rep(0, 200), bicameral = rep(1, 200),
                              federal = rep(0, 200), EUmember = rep(1, 200), growav = rep(2.455, 200),
                              xtimeyrs = rep(3.707, 200), opppreel = rep(0, 200), outalways = rep(0, 200),
                              npledgediv10 = rep(17.74, 200), precoagree = rep(0, 200), 
                              degomedlogrile = rep(0.2491240, 200), y1980s = rep(0, 200),
                              y1990s = rep(0, 200), y2000s = rep(1, 200), subset = rep(0, 200),
                              women_cabinet = seq(0, 0.6, length.out=200),
                              women_cabinet_sq = (seq(0, 0.6, length.out=200))^2,
                              partyid = rep(factor(701, levels(as.factor(df_tab2$partyid))), 200))
preds_wcabf1sq <- predict(m1bsq, newdata = newdat_wcabf1sq, se.fit = TRUE)
cis_wcabf1sq <- round(glm.cis(as.numeric(preds_wcabf1sq), sqrt(attr(preds_wcabf1sq, "var")), 0.95, m1bsq$df.residual),4)
preds_wcabf2sq <- predict(m2bsq, newdata = newdat_wcabf2sq, se.fit = TRUE)
cis_wcabf2sq <- round(glm.cis(as.numeric(preds_wcabf2sq), sqrt(attr(preds_wcabf2sq, "var")), 0.95, m2bsq$df.residual),4)
preds_wcabf3sq <- predict(m3bsq, newdata = newdat_wcabf3sq, se.fit = TRUE)
cis_wcabf3sq <- round(glm.cis(as.numeric(preds_wcabf3sq), sqrt(attr(preds_wcabf3sq, "var")), 0.95, m3bsq$df.residual),4)

pdf("cabinet_linear.pdf", height=8,width=6)
plot(newdat_wcabf1$women_cabinet, cis_wcabf1[,4], type="n",xlab="",ylab="",  yaxt="n", xaxt="n",
     ylim = c(0.3,1))
lines(newdat_wcabf1$women_cabinet, cis_wcabf1[,4], lwd=2, lty=1)
lines(newdat_wcabf2$women_cabinet, cis_wcabf2[,4], lwd=2, lty=2)
lines(newdat_wcabf3$women_cabinet, cis_wcabf3[,4], lwd=2, lty=3)
rug(df_tab2$women_cabinet)
axis(1, tck=0.03, cex.axis=1.5, mgp=c(0.3, 0.3, 0), lty=1, lwd=0, lwd.ticks = 1)
axis(2, tck=0.03, cex.axis=1.5, mgp=c(0.3, 0.3, 0), lty=1, lwd=0, lwd.ticks = 1, las=2)
title(xlab = 'Share of Women in Cabinet',
      ylab="Probability of Pledge Fulfillment",
      line = 2.25, cex.lab=1.5)
dev.off()

pdf("cabinet_sq.pdf", height=8,width=6)
plot(newdat_wcabf1$women_cabinet, cis_wcabf1[,4], type="n",xlab="",ylab="",  yaxt="n", xaxt="n",
     ylim = c(0.3,1))
lines(newdat_wcabf1sq$women_cabinet, cis_wcabf1sq[,4], lwd=2, lty=1)
lines(newdat_wcabf2sq$women_cabinet, cis_wcabf2sq[,4], lwd=2, lty=2)
lines(newdat_wcabf3sq$women_cabinet, cis_wcabf3sq[,4], lwd=2, lty=3)
rug(df_tab2$women_cabinet)
axis(1, tck=0.03, cex.axis=1.5, mgp=c(0.3, 0.3, 0), lty=1, lwd=0, lwd.ticks = 1)
axis(2, tck=0.03, cex.axis=1.5, mgp=c(0.3, 0.3, 0), lty=1, lwd=0, lwd.ticks = 1, las=2)
title(xlab = 'Share of Women in Cabinet',
      ylab="Probability of Pledge Fulfillment",
      line = 2.25, cex.lab=1.5)
dev.off()

#####
